import numpy as np
import pandas as pd
import plotly.express as px
import math
import plotly.io as pio
import plotly.graph_objects as go
# pio.renderers.default='notebook'
pd.options.plotting.backend = "plotly"
# utils
# def plot_state_vs_time(df):
# fig = df.plot()
# fig.update_layout(height=700)
# fig.show()
def plot_state_vs_time(df):
fig = go.Figure()
fig.add_trace(go.Scatter(
x=df.index,
y=df.iloc[:, 0],
legendgroup="group", # this can be any string, not just "group"
legendgrouptitle_text=df.columns.get_level_values(0)[2],
name=df.columns.get_level_values(1)[2],
mode="lines",
marker=dict(color="blue", size=10)
))
fig.add_trace(go.Scatter(
x=df.index,
y=df.iloc[:, 1],
legendgroup="group",
name=df.columns.get_level_values(1)[1],
mode="lines",
line=dict(color="red")
))
fig.add_trace(go.Scatter(
x=df.index,
y=df.iloc[:, 2],
legendgroup="group2", # this can be any string, not just "group"
legendgrouptitle_text=df.columns.get_level_values(0)[1],
name=df.columns.get_level_values(1)[2],
mode="lines",
marker=dict(color="green", size=10)
))
fig.add_trace(go.Scatter(
x=df.index,
y=df.iloc[:, 3],
legendgroup="group2",
name=df.columns.get_level_values(1)[1],
mode="lines",
line=dict(color="purple")
))
fig.update_layout(height=700)
fig.show()
$\underline{X(t_{k + 1})} = \underline{X(t_{k})} + dt\underline{\dot{X(t_{k})}}$
def euler_integrator(x_0, t_0, t_n, dt, dx):
x = x_0
t = t_0
yield (x, t)
while t <= t_n:
# eval x_k+1 using x_k and t_k
x = x + dt * dx(x, t)
# t increament to t_k+1
t = t + dt
yield (x, t)
https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
$\underline{X(t_{k + 1})} = \underline{X(t_{k})} + dt(k_1 + 2k_2 + 2k_3 + k_4)\frac{1}{6}$
where:
$f(x,t) = \underline{\dot{X(t_{k})}}$
$k_1 = f(\underline{x_k}, t_k)$
$k_2 = f(\underline{x_k} + dt\frac{k_1}{2}, t_k + \frac{dt}{2})$
$k_3 = f(\underline{x_k} + dt\frac{k_2}{2}, t_k + \frac{dt}{2})$
$k_4 = f(\underline{x_k} + dt k_3, t_k + dt)$
def rk4_integrator(x_0, t_0, t_n, dt, dx):
x = x_0
t = t_0
yield (x, t)
while t <= t_n:
# eval x_k+1 using x_k and t_k
k1 = dx(x, t)
k2 = dx(x + dt * k1 / 2, t + dt / 2)
k3 = dx(x + dt * k2 / 2, t + dt / 2)
k4 = dx(x + dt * k3, t + dt)
x = x + dt * (k1 + 2*k2 + 2*k3 + k4) * (1/6)
# t increament to t_k+1
t = t + dt
yield (x, t)
$\ddot{x}(t) = -2x(t) - 10\dot{x}(t)$
Let
$x(t) = x_1(t)$
$\dot{x}(t) = x_2(t)$
We define the state vector $\textbf{X} = [x_1, x_2]^T$
We can rewrite the equations given in the question as the following system of 2 first order O.D.E.s.
$\dot{x_1}(t) = x_2(t)$
$\dot{x_2}(t) = -2x_1(t) - 10x_2(t)$
A = np.array(
[
[0, 1],
[-2, -10]
]
)
def dx(X, t = None):
return A@X
X_0 = np.array(
[
[1],
[0]
]
)
T = 20
# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")
# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")
state_time_df = pd.concat([r_state_time_df, e_state_time_df])
plot_state_vs_time(state_time_df)
$2\ddot{x}(t) - 5\dot{x}(t) + x(t) = 0$
Let
$x(t) = x_1(t)$
$\dot{x}(t) = x_2(t)$
We define the state vector $\textbf{X} = [x_1, x_2]^T$
We can rewrite the equations given in the question as the following system of 2 first order O.D.E.s.
$\dot{x_1}(t) = x_2(t)$
$\dot{x_2}(t) = \frac{5x_2(t) - x_1(t)}{2}$
A = np.array(
[
[0, 1],
[-1/2, 5/2]
]
)
def dx(X, t = None):
return A@X
X_0 = np.array(
[
[1],
[0]
]
)
T = 20
# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")
# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")
state_time_df = pd.concat([r_state_time_df, e_state_time_df])
plot_state_vs_time(state_time_df)
$\ddot{x}(t) = 2y(t) - sin(4t^2x(t)), x(0) = 1, \dot{x}(0) = 2$
$\ddot{y}(t) = -2x(t) - \frac{1}{2t^2x(t)^2 + 3}, y(0)=1, \dot{y}(0) = 1$
Let
$x(t) = x_1(t)$
$\dot{x}(t) = x_2(t)$
$y(t) = y_1(t)$
$\dot{y}(t) = y_2(t)$
We define the state vector $\textbf{X} = [x_1, x_2, y_1, y_2]^T$
We can rewrite the equations given in the question as the following:
$\dot{\textbf{X}}(t) = f(\textbf{X}) = \begin{bmatrix}x_2 \\ 2y_1 - sin(4t^2x_1) \\ y_2 \\ -2x_1 - \frac{1}{2t^2x_2^2 + 3 } \end{bmatrix}$
def dx(X, t):
x_1 = X[0][0]
x_2 = X[1][0]
y_1 = X[2][0]
y_2 = X[3][0]
return np.array([
[x_2],
[2 * y_1 - np.sin(4 * t**2 * x_1)],
[y_2],
[-2 * x_1 - 1 / (2 * t**2 * x_2 ** 2 + 3)]
])
X_0 = np.array(
[
[1],
[2],
[1],
[0]
]
)
T = 2
# Euler
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = euler_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("Euler", "X_1"): x[0][0], ("Euler", "X_2"): x[1][0], "Time": t})
e_state_time_df = pd.DataFrame(state_time).set_index("Time")
# RK4
dt = 1e-3 # delta time (seconds)
state_time = []
state_integrator = rk4_integrator(X_0, 0, T, dt, dx)
for x, t in state_integrator:
state_time.append({("RK4", "X_1"): x[0][0], ("RK4", "X_2"): x[1][0], "Time": t})
r_state_time_df = pd.DataFrame(state_time).set_index("Time")
state_time_df = pd.concat([r_state_time_df, e_state_time_df])
plot_state_vs_time(state_time_df)